USHER: an algorithm for particle insertion in dense fluids 
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The insertion of solvent particles in molecular dynamics simulations of complex fluids is required 
in many situations involving open systems, but this challenging task has been scarcely explored 
in the literature. We propose a simple and fast algorithm (usher) that inserts the new solvent 
particles at locations where the potential energy has the desired prespecified value. For instance, 
this value may be set equal to the system's excess energy per particle, in such way that the inserted 
particles are energetically indistinguishable from the other particles present. During the search for 
the insertion site, the usher algorithm uses a steepest descent iterator with a displacement whose 
magnitude is adapted to the local features of the energy landscape. The only adjustable parameter 
in the algorithm is the maximum displacement and we show that its optimal value can be extracted 
from an analysis of the structure of the potential energy landscape. We present insertion tests 
in periodic and non-periodic systems filled with a Lennard- Jones fluid whose density ranges from 
moderate values to high values. 
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Many dynamical processes of chemical, biological and physical interest occur in open systems where matter and 
energy are exchanged with the surroundings. The main focus of attention has been in (Grand Canonical) Monte 
Carlo algorithms, which are particularly suited for the study of equilibrium states. However more recently there has 
£j ' been significant attention focused on molecular dynamics (MD) algorithms adapted for open systems. One of the 
biggest challenges faced in the investigation of such phenomena by MD simulation is the problem of efficient insertion 
of solvent particles in dense liquids. Indeed, as the scope and scale of MD increases, agrowing variety of applications 
and methods are in need of on the fly particle insertion algorithms (see Refs. 0, la, La 13 an d references therein). 
Among those, a particularly relevant family of methods are hybrid schemes that couple the particle domain to an 
\Q ' outer region described by continuum fluid dynamics 0, @ ■ The present method was devised for use in such hybrid 
C*~) t schemes, but we believe its application within molecular simulation may prove to be more widespread. 

The problem of inserting a solvent molecule in a dense fluid is commonly encountered in Grand Canonical Monte 
Carlo methods, for instance in Gibbs ensemble calculations for phase equilibria or evaluation of the chemical potential. 
A number of techniques have been proposed to overcome this problem (see 0,0 and references therein). For instance, 
cavity-biased procedures search for domains within the fluid with a small local value of number density, as these cavities 
are more susceptible to accommodate a new molecule. In doing so, a bias is introduced and, according to the rules 
of MC simulation, the bias has to be precalculated and corrected so that the scheme adheres to detail balance. Even 
so, recent comparisons showed that if the molecules are smaller than the mean size of the cavities, GCMC is nearly 
ten times faster than GCMD j| ■ We believe that the insertion algorithm proposed here may be used to improve the 
Ch \ efficiency of the GCMD schemes. 

The acme of a particle insertion protocol for MD is one that, in just a few iterations, is able to place the new 
particle within the required subdomain of the simulation space at a site where the potential energy takes exactly the 
desired value. This last condition ensures that no extra energy is introduced into the system, and therefore such an 
insertion algorithm would not require thermostatting after each insertion. Indeed, even for moderate liquid densities, 
these are difficult requirements and the few insertion protocols proposed in the literature 0, |(| are far from fulfilling 
' them. 

For instance, Goodfellow et al. 6] introduce solvent (water) molecules in the cavities of proteins to investigate their 
structural stability. Once the protein's cavities are found, the insertion protocol consists of several steps that involve 
operations over the whole system. Solvent molecules are introduced with arbitrary orientation and locations within 
the selected cavity. As a consequence, the energy of the system increases sharply after each solvent molecule insertion 
and, to allow its relaxation, two hundred energy minimization steps of the whole system (protein + water) are then 
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performed, followed by a one picosecond molecular dynamics simulation. This expensive insertion procedure, which 
involves substantial alterations of the microscopic dynamics, could be avoided if the solvent particle were initially 
introduced at the desired potential energy site. 

The work by Pettitt and coworkers, towards the Grand Canonical Molecular Dynamics method £| (see also [f| 
for further developent and applications), is an example of an open system MD simulation which does takes care of 
the potential energy at the insertion site. In their method, new particles contribute a fractional number to the total 
number of particles. These fractional or scaled particles must be inserted at positions where the potential energy is 
equal to that of the former (added or deleted) fractional particle. As explained in Ref. Q, the authors first use a grid 
method to slice the MD domain into a number of boxes that is the same as or a little larger than the total number of 
particles. The most favourable boxes (with the least number of neighbours) are selected as candidates to add the new 
particle. Then the new solvent candidates are placed within each of these boxes and two hundred possible molecule 
orientations of these new solvent candidates are computed. For each box, the orientation that yields the potential 
energy closest to the desired value is chosen and a first steepest descent procedure of O(10) steps follows. If this does 
not lead to any site with the desired insertion energy, they finally perform a large steepest descent procedure with at 
least 100 steps on the most favourable box. As the authors mention, this numerically expensive protocol still yields 
numerical errors that can disturb the system |4| . 

Our main concern is the insertion of solvent particles in the framework of a hybrid (particle-continuum) scheme. In 
recent work we proposed a hybrid shcmc that is able to deal not only with momentum but also with mass and energy 
exchange between the continuum (C) and the particle regions (P) [8j. In particular, particles need to be inserted in 
the overlapping C — > P regions where the C-fluxes are imposed on the P-domain. In a real liquid (with interacting 
potential energy), mass and energy exchanges are strongly coupled and we showed that, in order to balance the energy 
flux, the new particles have to be inserted at positions where the potential energy equals the value prescribed by the 
continuum domain. In that work yl we used a particle insertion algorithm (called usher) which is able to tackle this 
task in a rather efficient way (see 8] for a brief description). Further research has led to an enhanced version of the 
usher algorithm. Here we shall describe this new version of the USHER protocol and, for the sake of consistency, we 
shall briefly review the one one presented in Ref. 

The rest of the article proceeds as follows. We first formulate the root-finding problem in Sec. [H] In Scc lIIII we 
describe a reference scheme against which we compare the usher algorithm and then the USHER protocol. Insertion 
tests are described in Sec. II VI and the results are discussed in Sec. [3 In Sec. I VII we present an analysis of the 
potential energy landscape that proves to be very useful for the optimisation of the algorithm's parameters. Finally, 
conclusions and directions for future research on insertion algorithms are given in Sec. IVIII 

II. THE INSERTION PROBLEM 

We consider a set of N particles inside a box of volume V which interacts via pair potentials, V(r). At any instant 
the potential energy can be defined at any point r by evaluating U(r) = | J2iLo^(\ r — r i|)- The force that a test 
particle would feel at any point can be measured by f(r) = — VZ7(r). In this work, we consider a Lennard-Jones fluid 
whose interparticle potential V(r) = 4(r~ 12 — r~ 6 ) is written in the usual units of length (the effective radius a) and 
energy (the potential well e). 

The objective of the algorithms presented below is to find a position ro for which the potential energy equals a 
prescribed value, Uq] therefore U(r ) = Uq. In most practical situations a less stringent requirement needs to be 
fulfilled, namely (U(vq)) — Uq, where brackets denote an average over a certain (small) number of insertions. 

Even for a simple system such as the Lennard-Jones fluid, the structure of the energy landscape is very complex, 
with large energy gradients, and complicated energy isosurface shapes. A typical energy distribution along the whole 
space spreads over several order of magnitudes, but for the typical (moderate) temperatures usually considered in 
applications, the particles need to be placed at positions with extremely low energies compared with the range of 
the energy distribution. The result is that the mean specific excess energy resulting from the equation of state 
u eos — U eos (p, T) /N is a very low energy compared with the typical energies found at any arbitrary point of the 
space. As the fluid particle density increases the situation worsens, reflecting the fact that particles tend to reside 
within deep potential wells. The relation of the chemical potential p with density p shows clear evidence of this fact. 
At moderate densities the value of p is close to u eos but, above a certain density, p steeply increases above u eos , 
meaning that the typical energy needed to insert a particle becomes much larger than the mean potential energy per 
particle. 

Therefore if one needs to insert particles at positions with energies close to the mean excess energy per particle, 
u eos , one needs to find extremely low energy sites, particularly in dense systems. The main problem to be faced 
is that the energy landscape presents many energy " holes" whose local minima range from intermediate to large- 
energy values. Here, we define a "hole" as a region of space enclosed by an isosurface of energy in such way that 
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VJ7(r) • n > at the hole surface, where n is the (outward) normal surface vector. Usually, these holes act as traps 
for the widely used energy-minimisation algorithms based on the standard steepest descent or conjugate gradient 
methods 0- As a matter of fact, we soon discovered that it was very inefficient to move downhill over the potential 
energy by means of any of the standard versions of the steepest descent method used in molecular simulations (see 
for instance, Ref.0). The purpose of the present study is to present a stcepest-descent-like iterative procedure that 
can sort out the intermediate-energy holes. To this end the USHER algorithm does not rely on line minimisation along 
the steepest descent direction [Tol |. but instead on a displacement size which is adapted on the fly, according to the 
local topology of the potential energy landscape. Another advantage is the facile implementation of the usher code. 

III. DESCRIPTION OF THE ALGORITHMS 

We shall now describe some general aspects of the problem of particle insertion and the common features of the 
algorithms concerned. 

In any insertion procedure the first step is to place the new particle at a starting position r' '. In all the tests 
presented in Sec. IIVI we chose r' ' at random. We also tried to select r^ ) according to a cavity-biased procedure 
(as in Ref. 3). As explained in Sec. this procedure incurs a number of operations of O(N) prior to the insertion 
algorithm itself and we found that, when using the usher algorithm presented below, it did not reduce the total 
number of iterations with respect to the (much cheaper) random choice. 

During successive iterations the iterator's position is moved according to the update rule which, in general, may 
be a function of the mechanical quantities at the previous iteration, r( n+1 ) = r( n+1 ) (r^ n \ U^ n \ f( n )) . The search 
terminates if the new position r(" +1 ) is a site with the desired energy. This is determined by the following condition 

l£l ( " +1) <£ max , with = u{n Z 7 Ua (1) 

\ u o\ 

where £ max is a pre-determined parameter, namely the half-width of the interval of the accepted energies around Uo, 
and £ is defined as the relative difference of the potential energy at the (n + 1) iteration with respect to the 

desired value Uq. 

Finally, once the new particle is correctly inserted, the force that it exerts on its neighbours is calculated and its 
velocity is also assigned. This velocity is drawn from a Maxwellian distribution with the desired temperature T and 
the desired mean velocity, (v) 

P W=(wr) CXP ( 2mkT < 2 > 

While the algorithm is guiding a new particle to a correct location, the positions of all the other particles remain 
frozen. This means that one insertion iteration only involves the evaluation of the force on a single particle (that is, 
the force exerted by all the particles at the site r^'^). 

The starting position determines whether the following iterations will have to be downhill (if > Uo) or uphill (if 
< A simple way to unify both cases in a single scheme is to rescale the potential energy as U (r) — > sgn U(r), 

where sgn = ^| )Z^° ■ By doing so the forces f = — Vf7 are also redefined and, in particular, a case with sgn = — 1 

then implies that the redefined forces point uphill of the (unsealed) potential energy throughout the entire course of 
that particular particle insertion. In the following presentation we shall assume that the energy and force field are 
already rescaled as sgnU(r) and sgnf, so we will not explicitly include sgn in the equations. 

Durin the iterative process the algorithm will encounter three different situations which may require a separate 
treatment (for instance, different update rules). We denote these situations as follows as: downhill move, < 
UW; uphill move, [/( n+1 ) > {/(") and confinement, [7 < -™ +1 ^ < Uo < U^ n \ In an optimal insertion one expects to keep 
going downhill (with respect to the rescaled potential energy) until the confinement is attained. Then, the desired 
location lies within the segment Sr^ = r^" +1 ) — rW and can be determined by means of standard one-dimensional 
(ID) root-finding algorithms (such as the Newton-Raphson or bisection methods). 

The most problematic iteration corresponds to the uphill move and it deserves some discussion. To illustrate the 
scenario we refer to the energy landscape shown in Fig. 1. Even for moderate densities (Fig. 1 
corresponds to p = 0.6), the low energy regions conform to a complex tube-like structure. The insertion algorithm 
will have to usher the new particle into these energy tubes before arriving at a correct location. An uphill iteration 
may arise when the iterator faces either of two features of the energy landscape: intermediate-energy holes or sharp 
bends (including saddle points). Note that both kind of features induce completely different decisions. The best thing 
to do when encountering an energy trap is to give up the search and restart from another initial position r^ ' . By 
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contrast, if a bend in the energy landscape leads to a low energy valley, it may be worthwhile to use an update rule 
that can efficiently deflect the iterator's trajectory. Unfortunately, once an uphill move occurs it is not possible to 
distinguish between these two features within only one iteration. On the other hand, the number of uphill iterations 
rapidly increases as the displacement |r(™ +1 ) — r^| is made larger than a specified maximum threshold. In fact, an 
important issue for the algorithm design is first to estimate this threshold and then, to determine the best decision 
to take upon an uphill move (see Sec. IVIJI . 

It is also convenient to introduce a restart condition in order to avoid any possible stagnation of the algorithm 
around energy holes. In particular, if n > n max the search is restarted from another initial position r^ ). Particularly 
at high densities (typically above 0.75), the overall number of iterations is sensitive to n max . A very large value of n max 
corresponds to many unsuccessful and time-consuming iterations, while a value for n max that is too small prevents 
most of the potentially successful trials from terminating successfully. We found that the best compromise between 
these two extremes is to make n max ~ 0.8 (n), where (n) is the number of iterations averaged over a certain number 
(~ 20) of insertions. The value of n max depends on the density. It may be determined from an initial test-run, or 
alternatively reassigned on the fly according to the value of (n) determined during the simulation. 

In the remainder of this section we first define a "reference" scheme against which we can then discuss and compare 
the USHER algorithm in more detail. 

A. The reference scheme 

In order to better understand the behaviour of the usher algorithm it is helpful to compare its performance with 
a reference scheme based on a combination of well established methods widely used in the literature for root-finding 
and energy minimisation. While moving downhill, the reference scheme uses a basic steepest descent step with a fixed 
displacement As\. The update rule being is 

f (n) 

r (n+l) = r (n) + A*, (3) 

where, according to standard notation, / W is the modulus of f (™) . 

If an uphill move is made, the reference scheme will first try to deflect the iterator's trajectory in order to adapt 
itself to a possible bend in the potential energy surface. By construction of the update rule, Eq. (0, the potential 
energy decreases locally at r™' in the direction Sr^ = r( n+1 ) — r^ n \ Therefore if [/( n+1 ) > f/W there must exists a 
location r m where, U(r m ) = min {t/(rA)|rA = r( n ) + A5r'™H. The reference scheme finds the position r m by means 
of a line-minimisation of the potential energy along the segment 5r^ (see [T^j f° r details). The new position is then 
recalculated by a steepest descent step starting from r m and with a displacement As2; 

r(" +1 )=r m + ^A S2 . (4) 

The line-minimisation itself requires an inner iterative procedure (see [Tol|'). In view of the narrowness of the potential- 
energy tubes, we used no more than three iterations for the estimation of r m . Better estimates of r m do not improve 
the efficiency, but instead lead to a larger number of force evaluations in the overall scheme. When a local minimum 
of potential energy is found, the iterator's position will bounce back, moving subsequently upwards and downwards in 
energy as it hits the walls of the energy hole. To avoid this situation, after several (typically three) consecutive uphill 
iterations the scheme determines that it has been trapped in a local energy minimum and consequently restarts the 
search from another initial position r^ ^. Finally, once the root has been confined in the segment Sr^ n ' — r( n+1 ) — r^ n ', 
it is expressed as ro = r( n+1 ) + XqSt^, with Ao a real number in (0,1). To find Ao, the reference scheme uses a 
ID root-finding algorithm which combines the Newton-Raphson method with the robust bisection method to ensure 
confinement in case of a failure of the Newton-Raphson step (due, for instance, to /w ~ 0). It typically took less 
than three iterations to calculate the value of ro. The optimum choice for the parameters Asi and As 2 is presented 
in Sec. \VJ\ 

B. The usher scheme 

The basic idea of the usher insertion algorithm is to use an update rule to move downhill that can adapt the 
iterator's displacement according to the local topology of the low energy landscape. This is reflected in the following 
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update rule: 



f (n) 

r (n+l) („) + J_5 S («). ( 5 ) 



Equation JSJ is essentially a steepest descent scheme with an displacement Ss^ that depends on the iterator's position. 
The success of the method resides in a judicious choice of 8s^ n > . Optimal performance was obtained using the following 
expression for Ss^ n \ which depends on both the local potential energy f/w and force fW, 

A if f/W > Uovlp 

^ {n) = < „:_/V UW-U \ :CTT(n) _ r (6) 



min (As, (n) ° j , if UW < U ovlp . 



The best way to illustrate how the adaptive displacement of Eq. JfjJ) works is to describe how the USHER scheme 
performs one insertion. As long as the starting position r^ ^ is chosen at random, there is a large chance of overlap 
with a pre-existing particle, leading to a very large value of U^°\ The displacement As vip quoted in the first line 
of Eq.JfJJ can be constructed to remove the overlap in (typically) one iteration. For this reason, U ov i p is chosen to 
be a very large energy representing an overlap position, say U ov i p ~ 10 4 . As the hard-core part of the interparticle 
potential goes like 4r~ 12 , the distance from a site with energy f/W > U ov i p to the centre of the overlapped particle 
is r = (4/C/(™)) 1 / 12 . Therefore by choosing As ov i p — r a — (4/C/W) 1 / 12 , we can guarantee that the next iterator's 
position r(™ +1 ) will be moved a distance r a away from the centre of the overlapped particle and, by virtue of Eq. 
in a direction of lower potential energy. The value of r a should be close to or slightly smaller than a characteristic 
contact distance between particles (e.g. the distance given by the maximum in the radial distribution function). For 
the pure Lcnnard- Jones fluid under consideration here, we have used r a = 0.9 (in units of a). 

Once any possible initial overlap is sorted out (U < U ov i P ), the second line of Eq. © is designed to drive the 
new particle downhill in energy, towards the target value Uq. Here resides the main difference with respect to the 
reference scheme. At large energies, the typical slope of the potential energy is very large (/w >> 1), meaning that 
the energy drop along the steepest descent direction is governed by the linear term of the Taylor expansion in the 
displacement, AU = f^ n 'Ss + 0(5s 2 ). The second line of Eq. J§J makes use of this fact and takes AU = f/W — Uo for 
extracting a displacement adapted to the (maximum) local energy gradient Ss = AU //W. Note that at large energies 
[/(«•) _ u Q f/W, so after one iteration one expects the energy to decrease in (at least) a fraction of U^ n ', meaning 
a linear convergence. The local curvature of the potential energy landscape becomes dominant when approaching a 
local minimum (/W - 0) and in this case Eq. © limits the displacement to a maximum value As. The maximum 
displacement is the only variable parameter in the algorithm and, as discussed in Sec. I VII its optimal value is about 
the width of the low-energy tubes of the potential energy landscape (see Fig. 1) 

At low energies, as the iterator approaches the energy target Uq, the displacement Ss — (f/W — Uo)/f^ behaves 
like a Newton-Raphson step made along the steepest descent direction. Due to this feature, the convergence of the 
usher algorithm increases notably near the target. In particular, this kind of displacement enables the error £ to 
decrease quadratically once £ < 0(1). This fact is illustrated in Fig. 2 by plotting the absolute value of the error 
£|(n+i) a g a i ns t its value at the previous iteration |£|W. As explained above, for £ > 0(1) the algorithm converges 

linearly with |£|("+ 1 ) ~ 0.4|fW| ; while |£| (n+1) ~ 0.35 (£ (n) ) 2 once £ < 0(1). In the same way, it may be possible 
to further increase the con verg ence rate by implementing a displacement based upon higher order methods such as 
Halley's or Bailey's scheme [Tl| (for such purpose one would need to calculate the Hessian matrix and project it onto 
the steepest descent direction). 

For the sake of completeness, we also describe here the older version of the usher's displacement used in Ref. 
This earlier version used a similar displacement rule for f/W > U ov i p to that quoted in Eq. but for lower energies 
it used Ss^ — min(As, ^f^At) where the optimal choice for the parameter At ranged within (0.05,0.15). This 
scheme is around two times slower than the improved USHER sheme discussed here. 

One of the important issues of the algorithm design concerned the optimal strategy to deploy for uphill iterations. 
We compared two different strategies. The first one, which we shall call indirect usher performs a line-minimisation 
[loj of the energy along the direction Sr^ = r(' i+1 ) — rW, similar to that described in Sec. IIII Al and Eq. The 
second alternative, called direct-VSHER, gives up the initial search and restarts a new one from another random position 
r(°) once an uphill move is encountered. Interestingly the insertion tests (see Sec. IIV|I clearly show that the direct 
usher is about two times faster than the indirect version. This indicates that most of the uphill moves encountered 
using the update rule of Eqs. (JSJ and JBJ) are due to energy holes and therefore suggests that Eq. (JBJ enables the 
usher algorithm to properly deflect its trajectory at most of the bends of the low energy-tubes encountered. A less 
restrictive version of the direct usher allows a line-minimisation iteration only if the uphill move is done near enough 
to the target (for instance if |£| < O(l)). This alternative gives slightly better results at large densities. 
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In the insertion tests presented below in Sec. IIVI the reference scheme is compared with the most efficient version 
of the USHER algorithm; i.e. with the direct usher. To avoid any possible confusion, in the remainder of the paper 
we shall simply called this, the usher algorithm. 

IV. INSERTION TESTS 

The insertion algorithms presented in Sec. IIIII were evaluated in two kinds of systems, with and without periodic 
boundary conditions. We stress that no thermostat was used in any of the insertion tests. This ensures that the 
temperature of the system does not spuriously increase due to the dissipation of possible additional internal energy 
introduced by particle insertions in non-appropriate (higher-energy) locations. 

In order to investigate the functioning of the insertion algorithm we shall drive the system through a specific 
thermodynamic process (see below) and compare the values of the thermodynamic variables computed during the 
simulations with those arising from thermodynamics. The system contains N particles within a volume V and its 
total energy is E = 3NT/2 + U(r N ), the energy per particle being e = E/N. The thermodynamic processes will be 
specified by the variation of the number of particles AN (or density Ap) and the change of energy per particle Ae. 
We now use standard thermodynamics to derive the changes in the system's other variables. 

The variation of energy per particle upon insertion of AiV particles into the system is Ae = AE/N — eAN/N. 
For a system having no contact with a thermostat or a manostat, as for the one considered here, the variation of the 
total energy upon insertion of AiV particles is exactly AE — (e 1 ) AN, where (e') is the energy of the inserted particle 
averaged over AN insertions; AiV (e') = ( e', being the energy of the ith inserted particle). Thus, 

Ae = « e '> - e) AN/N = «e') - e) Ap/p. (7) 

The variation of temperature can now be calculated from the equation 

e -(l),H (8) 

where Ap — (AN)/V, c v is the specific heat at constant volume and (de/dp)r is obtained from the equation of state 
for the excess energy per particle u = U/N reported by Johnson et al. 

Therefore, for given initial values of the system's density and temperature (pq,Tq), the time evolution of the 
thermodynamic variables is determined by the (specified) temporal variation of density dp/dt. The rate of temperature 
variation, obtained from Eq. ©, enables us to calculate the temperature at each instant in the process. The pressure 
and the excess energy per particle can be then obtained from equations of state (P — P eos (p, T) and u — u eos (p, T)). 

In the tests presented below we considered a thermodynamic process in which the density increases at constant 
specific energy Ae = 0. According to Eq. Q, during the process the average energy of the inserted particles (e') 
is set equal to the mean specific energy of the system, e. This condition is similar to that required for the energy 
balance conditions in the hybrid (particle-continuum) scheme of Delgado-Buscalioni and Coveney Q. In fact, the 
process with Ae = can also be sought as a test for energy conservation in this hybrid scheme. 

The thermodynamic relations, such as Eq. ©, are meaningful at least under condition of local equilibrium. This 
imposes a limit on the rate of particle insertion, because within each subdomain of the system the insertions of particles 
need to be sufficiently well spaced out in time for the system to be able to recover the equilibrium distribution. Consider 
a small subvolume of size A, large enough to be representative of the system's distribution function. For instance A 
may be the distance at which the radial distribution function converges to one (~ 3<r). To be sure that the system is 
able to restore its equilibrium distribution, the rate of particle insertion in each of these subdomains A 3 (dp/dt) has 
to be smaller than the inverse of the collision time, which for a simple fluid can be estimated using the hard-sphere 
approximation by t c ~ 0.14p~ 1 T~ 1 / 2 (cr 2 m/e) 1 / 2 (see e.g., 0). In our calculations we used (dp/dt) < 0.01, so the 
characteristic insertion time was ~ 3, much larger than the collision time, t c ~ 0.3 

A. Insertions in a periodic box 

The first set of tests were performed in systems contained within a cubic periodic box of side lengths L — {7, 8, 10} a. 
The initial density was set to a moderate value p(t — 0) = 0.4 and was increased until p ~ 1.0. The maximum rate 
of density increase used was dp/dt ~ 0.01. The temperature, pressure, excess energy per particle u — U/N and total 
energy per particle e are plotted in Fig. 3 versus the density. Results correspond to particle insertions in a box with 
L = 10a at a constant density increase rate of dp/dt = 0.01. Particles were inserted at sites where the potential 
energy equals the specific excess energy of the system Uq = U /N = u and with velocities drawn from a Maxwellian 
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distribution at the instantaneous system's (kinetic) temperature T. The dashed lines in Fig. 3 correspond to the 
thermodynamic variables obtained from the the equation of state according to the process of Eq. iJSJ . 

B. Insertions in an open flow 

The second test was done in open fluid flows, i.e., in systems with open boundary conditions. We considered a 
system with density N/V within a cuboidal domain of sides L x — 40 a and L y = L z = 9 a. The system is periodic in 
the y and z directions and has open boundaries at x = and x = L x . Particles were inserted with potential energies 
close the specific excess energy U eos /N and with velocities drawn from a Maxwellian distribution at a temperature 
To which was fixed throughout the simulation. Insertions were done within a region of width Ax around x = and 
at a rate Apvi n , where Vi n is a parameter that determines the flow velocity normal to the surface vector of the open 
boundary (n) and A — L y L z is the area of the boundary. At the right hand boundary, particles are extracted at the 
same rate, so the overall density of the system remains constant throughout the simulation N/V. In order to couple 
the particle region to the outer pressure we used our hybrid particle-continuum scheme at the x = and x = L x 
surfaces @- 

Note that in this case particle insertion in the x-direction is restricted to a region Ax which is set to Ax = 2. Oct (the 
volume available to insert particles being AxL y L z ). To ensure that insertions are done in this region two different 
strategies were implemented. The first one is a simple reflection of the position r(™ +1 ) back to the insertion domain 
when the usher iterator crosses the x — and x — Ax boundaries. In the second implementation we imposed an 
artificial sharp potential well at x = and x = Ax which acts only during the evaluation of the forces in the iteration 
procedure, i.e., it was not included in the evaluation of U(r^ n+1 ^). Both alternatives worked equally well and resulted 
in a similar number of required iterations. 

The simulation starts from an initial state with zero mean velocity and constant density profile along the x-direction. 
As time goes by, the particle insertions concentrated in the region around x — lead to the production of a density 
wave that expands at the sound velocity for x > 0. This density wave transports momentum along the x-direction and, 
after a transient time, the density profile converges to the flat stationary density profile; throughout the simulation 
cell, the mean flow x- velocity tends to the value vi n . 

The hydrodynamic and thermodynamic variables were measured over slices of width Ax along the x direction. 
Figure 4 shows the local density at some of the left most slices x < L x /2 together with the mean (slice averaged) 
velocity and the total temperature of the system. The oscillatory behaviour of the local density is a desired feature 
of these tests as it enables us to determine the dependence of the number of iterations n on the density for a range 
of values of p in each simulation. We refer to our previous paper 8] for a detailed comparison between theoretical 
hydrodynamic trends and results obtained from hybrid continuum-particle simulations in different relaxing flows, also 
involving mass exchange. 

V. RESULTS 

Figure 5 presents the average number of single- force evaluations (rif) needed to perform an insertion versus the 
density. The most expensive part of the insertion algorithms is the evaluation of the force. Consequently in order to 
compare the performance of the algorithms we have used rt/, rather than the total number of iterations n. We note 
that some of the steps discussed in Sec. IIIII fsuch as line-minimisation) require several sub-iterations and so n < n/. 
The USHER and the reference algorithms are compared in Fig. 5, in a test corresponding to insertions in a periodic 
box. In this kind of test (and for a similar average error (£) < 0.05) the usher algorithm is more than two times 
faster than the reference scheme for p > 0.5 and more than four times faster for p > 0.8. The reference scheme is 
slightly slower when insertions are constrained to a smaller region, as occurs in the open fluid flow tests. But notably, 
for both open flow and periodic box tests, the USHER scheme gives similar values of (nf) . This means that the usher 
algorithm does not pay any extra cost for restricting the size of the domain of insertion. This may be understood by 
looking at the distance between the initial trial and final insertion positions, Ar — |r(°) — r(™)|, shown in Fig. 6. For 
a wide range of densities the maximum value of Ar is smaller than l.Ocr (its average being typically less than 0.5er), 
indicating that most particles are inserted before the USHER iterator reaches the boundaries of the insertion domain. 
This feature important for many applications. For instance, in the hybrid particle-continuum schemes the insertions 
are assigned (and restricted) to finite cells arising from a discretisation of the space [?, 8] ; and in the water- insertion 
method [6|, the water molecules have to be placed in assigned protein cavities. 

Figure 5 illustrates how the number of force evaluations varies with the maximum averaged error when using the 
usher algorithm. In particular we compare the results for the insertions tests done at an initial temperature of To = 3 
(as in Fig. 2) and with different values of the maximum error averaged over 30 insertions, (£). To decrease the error, 
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from 0.15 to 0.03 one typically needs one more iteration. Another iteration leads to (£) = 10~ 3 . This fast (quadratic) 
error reduction is made possible by the Newton-Raphson-like displacement implemented in Eq. © (see Fig. 2). A 
systematic non-vanishing value of (£) has a direct effect on the thermodynamic variables, as shown in Fig. 3. For 
instance, a value of (£) ~ 0.05 maintained during the insertion process leads to a systematic drift from the Ae = 
line, and also has an effect on the temperature evolution. 

Additionally, Fig. 5 illustrates how (n/) varies with the system's temperature as shown in the data for To = 3 
and To = 10. At larger temperatures it becomes much easier to insert particles once p > 0.6. The reason is that the 
target energy Uo(= u eos {p,T)) increases much faster with the temperature at larger densities than it does at lower 
densities. For instance, u eos (0.4, 3) ~ —1.9 and u eos (0.4, 10) ~ —1.2, while for a larger density u eos (0.85,3) ~ —3.1 
and M eos (0.85, 10) ~ —0.6. For the same reason, if particles were inserted with potential energies similar to the 
chemical potential (Uq) — p eos (p,T), the slope of («/) with p would be flatter than those data shown in Fig. 5. We 
also performed insertions at subcritical temperatures T < 1.3, for liquid densities and also inside the liquid- vapour 
coexistence region. In these calculations the number of iterations needed to insert a particle was very similar to that 
presented in Fig. 5 for To = 3. Anyhow, fluctuations of n/ were larger inside the coexistence region as a consecuence 
of the inhomogeneity of the density field. 



VI. CHOOSING AN OPTIMAL PARAMETER SET 



We now wish to provide a physical interpretation of the performance of the particle insertion algorithms, based 
on the structure of the energy landscape. Such insight will be very useful for extensions of the usher parameters 
to the simulation of other kinds of fluids. In fact, our experience is that, instead of a simple parametric study, it 
is advisable perform an analysis of the structure of the potential energy landscape to obtain information about the 
typical shapes and length scales of the low energy regions. This kind of structural analysis for the Lennard- Jones 
fluid considered here not only provided important clues for the algorithm design, but provided the key relationship 
between the optimal displacement As and the density. 



A. Low energy holes 



In order to investigate the structure of the low energy holes we devised the following procedure. In a standard MD 
simulation in a periodic box and at time intervals separated by several collision times, we seek a point ro with a very 
low pre-specified energy Uq. In particular Uq is chosen to be the mean excess energy per particle. Initially the search 
for the point ro was done by the "basic" update rule of the usher algorithm mentioned in Sec. IIII Bl Once ro was 
found, the energy landscape was probed in radial directions from this point. For each azimuthal angle ip £ [0, 2n] and 
longitudinal angle 6 £ [0, w], the energy U(r) was measured for increasing radial coordinate and a radial distance was 
recorded when U(r) > Ui SO - The radial distance will be denoted as R(ip, 6). Note that R(ip, 9) determines the shape 
of each hole; in particular the mean radius and the mean of the squared radius were computed for each hole: 

i R )e,4> = A I* *!> r R{QA)d9, (9) 
Z7r Jo Jo 

( R \^^£ T dd>£R 2 (0,4>)dB. (10) 
The effective shape of the hole can be estimated by the following quantity, 

Clearly for a sphere or — 0, while (Jr is positive for any other elongated shape. A glance at the low and intermediate 
energy regions of a typical contour plot of the potential energy (see Fig. 1) suggests that it is possible to estimate 
the characteristic length scales of the low energy regions by fitting ur and (R) s ^ to ellipsoids. In particular, due 
to the symmetry of the LJ fluid, it is enough to use asymmetric ellipsoids for this estimation. For an ellipsoid with 
semi- minor and semi-major axes given respectively by R s and Ri — e R s , the following parametric relations fit within 
1% to the exact analytical results: 

a fl = 0.561ogx, (12) 
(R) e , = R S (1 + 0.25 log X ). (13) 
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For given values of <jr and (R) g ^, one can estimate the eccentricity \ = Ri/Rs and the semi- minor axis R s using 
Eqs. l(T2*f and l(T3|l . The values of (R) g rf) , <tr and the estimations of R s and Ri, averaged over a set of about 80 
holes, are shown in Fig. 7 versus the density. To ensure that these values are representative of the shape of low and 
intermediate-energy regions (such as those shown in Fig. 1), the averages were obtained for a relatively wide range 
of (intermediate) energy isovalues Ui SO G [0,50]. The error bars determine the maximum variation of these quantities 
for this range of Ui SQ . 

In Fig. 7 we included the optimum choice for the reference algorithm described in Sec. IIII Al fAsi and As 2 ) ■ The 
optimisation of these parameters was performed independently, for a wide range of densities, p 6 [0.4 — 0.92] (see Fig. 
8). It is interesting to note that Asi closely follows the trend obtained for the smallest effective radius R s , while As 2 
lies above the longest radius Ri of the intermediate energy regions. The interpretation of the results for Asi seems 
quite evident, meaning that the displacement of the steepest descent method when moving downhill should be about 
half the minimum characteristic diameter of the low energy valleys. From this we readily understand why the best 
choice for the maximum displacement of the usher algorithm is As ~ R s ~ 0.1p~ 15 . Quite remarkably, the estimate 
of the mean free path based on the hard-sphere fluid, 0.2p _1 (shown in Fig. 7), is close to the typical radius of the 
low energy holes (R) s ^. This indicates that such kinetic information, if available, may be of great help for the first 
adjustment of the maximum displacement of the algorithm, when inserting particles (or minimising the energy) in 
other kinds of fluids. 

Fig. 8b sheds light on the interpretation of the result for As 2 . The optimal value of As2 may be taken to be any 
value larger than a certain threshold, which according to Fig. 8b has to exceed the largest typical longest diameter 
within the low energy regions. This confirms that once an uphill move is made the fastest option is to completely 
traverse the energy valley and continue the iterations from a high energy site, instead of trying to pursue possible 
further line minimisations. This conclusion, obtained from the reference scheme, suggested that the best procedure 
was to give-up the search once 

jj(n) > fj(n-i)_ As sta ted in Sec. ImBl this is indeed what we have found when 
comparing the direct and indirect versions of the USHER scheme. 

VII. CONCLUSIONS 

An increasing number of methods involving molecular dynamics (MD) simulations of open systems 0, H, El 
require one to insert particles at precise locations where the potential energy is set equal to a pre-specified value. 
Moreover, insertions need to be done on the fly and the performance of these methods will greatly depend on the 
efficiency of the insertion algorithm. At high densities this may seem a formidable task and indeed this sort of insertion 
algorithms has scarcely been explored in the literature. The main purpose of our paper is to show that this problem 
can be efficiently accomplished. To this end we have devised a particle insertion procedure called the usher algorithm. 
To give an example, to insert a particle in a Lennard- Jones fluid with p = 0.5 and T = 3.0, at positions where the 
potential energy equals the mean specific energy of the system, the algorithm requires around 8 extra evaluations of 
a (single-particle) force and 25 if p = 0.8. 

The usher algorithm essentially consists of a steepest descent iteration procedure (see Eq. [5J with a displacement 
adapted to the local shape of the energy landscape. In particular, by using an initial displacement which depends on 
the value of the potential energy at the initial trial position, 1 — (4/C/(°)) 1 / 12 , the algorithm avoids in (about) one 
iteration any possible overlap with a pre-existing particle. We confirmed that this feature makes it advantageous to 
choose the initial trial position at random, instead of using a much more expensive grid method to slice up the entire 
space in the search of the less dense region (as is done in the cavity-biased Monte-Carlo or in the Grand-Canonical MD 
scheme proposed by Pettitt and coworkers 0,0). In subsequent iterations the displacement is given by the Newton- 
Raphson step measured along the steepest descent direction and has a upper bound of As, to avoid uncontrolled 
jumps near local minima. 

As another relevant conclusion, we wish to caution about the usage of line minimisation, normally implemented in 
conjunction with the steepest descent method 0,0. We clearly observed that, in these complex landscapes, it is 
better to use a (small enough) maximum displacement to ensure that most the iterations are made downhill; then, 
if a single iteration is made uphill, the best option is to restart the search from another random position, rather 
than performing a line minimisation. There are two reasons for this fact: firstly, line minimisation is expensive and 
secondly, and more importantly, we observed that if the maximum displacement is optimal, most of the uphill moves 
are due to the presence of an energy trap (i.e., a local minimum at an energy larger than the target). 

An important part of this work was to give a physical interpretation of the optimal maximum displacement As. 
As stated earlier it is essential to optimise its value for the complex topology of the low energy-tubes. To that end 
we analysed the structure of the energy landscape of the Lennard-Jones fluid considered here and concluded that 
the width of low-energy tubes scales as Q.lp~~ . A much more computationally expensive parametric study clearly 
showed that the optimal displacement follows a similar trend. This means that, in order to extend the insertion 
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algorithm to other fluids, it is strongly advisable to first investigate the structure (in terms of shape and length 
scales) of the low energy regions. Apart from this, our insertion protocol is based uniquely on mechanical variables 
readily calculated in any standard MD simulation (force and potential energy), therefore the same kind of protocol 
can be used for inserting solvent molecules in fluids consisting of (small) poly-atomic molecules or even in polar fluids 
with nonadditive potentials. Nevertheless the optimal algorithm design may require modifications, depending on the 
specific molecule considered. For instance, in fluids whose molecules have rotational degrees of freedom, the insertion 
update step could be modified so as to first update the position of the center of mass of the molecule and then to use 
the local torque to orientate the molecule to the most favourable position. Such an investigation will form the subject 
of future work. 
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IX. FIGURES 

Figure 1 

(a) A cut along the x=0 plane of the contour plot of the potential energy landscape for energies lower than 100, 
showing the typical low-energy tube-like structures. In (b) a close up of the left most region indicating with thicker 
solid lines a possible targeted energy, at U = —4. Some bends and saddle points of the energy surface and some energy 
traps (local minima with energy larger than the target) are indicated with solid and dashed arrows, respectively. The 
snapshot corresponds to a LJ fluid with p = 0.6 and T = 2.5 inside a 3D cubic box of side L — 10 a. The size of the 
maximum displacement As used by the usher algorithm for this density is also indicated. 

Figure 2 

The absolute value of the error at the n + 1 plotted against its value at the previous iteration n. The data corresponds 
to insertions made by the usher algoritm in a periodic box L = 10 a and for 0.6 < p < 0.75. As illustrated by the 
dashed lines, the convergence is linear for £ > O(l) and quadratic for £ < O(l). 

Figure 3 

(a) Total energy per particle e, (b) pressure P , (c) temperature and, (d) excess energy per particle u versus the 
density, obtained in a particle insertion test made in a cubic periodic box with side length L — 10 a. The density 
increases linearly with time at a rate of dp/dt = 0.01. The insertions were made to guarantee a process with Ae = 
(see text). The dashed lines corresponds to the thermodynamic variables extracted from Eq. JSJ, using the equations 
of state for u eos (p,T) and P eos (p,T). 

Figure 4 

Evolution in time of various hydrodynamic variables in an insertion test on an open flow in a "box" with sides 
L x = 40(J, L y = L z =9cr. Particles were inserted with zero mean velocity at x = at a rate s = L z L y poVi n and 
extracted at x = L x at the same rate. The overall density was po = 0.5 and Vi„ = 1.0. (a) Mean and local temperature 
and (b) local density at slices of width Ax = 2; in (c) the local velocity at each slice (dotted line) and the mean 
velocity (solid line) are shown. The flat stationary :r-profiles of density (p — 0.5) and velocity (v x — Vi„ — 1.0) are 
reached after several sound transversal periods. 

Figure 5 

The average number of force evaluations needed to insert a new particle (n/). The results correspond to insertions in 
a cubic periodic box of side L — 10 a. For all the curves (£) indicates the maximum value of the averaged errors and 
the error bars corresponds to the standard deviation upon 100 insertions. The results are for processes with Ae = 0, 
starting from an initial temperature To [see Eq. @ and Fig. 3]. 

Figure 6 

The distance travelled by the usher algorithm between the initial trial position and the final insertion site, Ar = 
| r (°) _ r'"'|, (in log-scale) versus the density. The test corresponds to particle insertions in a cubic periodic box of 
side L = 10 a. 

Figure 7 

The mean radius of the low energy regions {R) g ± [Eq. along with the smallest R 3 and largest Ri characteristic 
lengths of the low energy regions estimated by Eqs. 11 ol and l|12lL Squares correspond to the optimum values of the 
reference scheme Asi obtained from a parametric study. The dashed line (O.lp ) corresponds to our choice for 
the optimum usher maximum displacement As. The mean free path (hard-spheres estimate 0.2P" 1 ) is also shown. 
The inset shows the normalised variance an given by Eq. il l H and the estimated mean eccentricity of the low-energy 
holes x= Ri/Rs- 
Figure 8 

The average number of force evaluations per insertion versus the parameters (a) Asi and (b) As2 of the reference 
scheme of Sec. IIII A1 The evaluations were done by inserting particles in a cubic periodic box of side L = 10 a. The 
range of densities at which the evaluations were made is indicated in each figure. 




y 



FIG. 1 . R. Delgado-Buscalioni and P . Coveney 




Fig2 . R. Delgado-Buscalioni and P. Coveney 




FIG. 3. R. Delgado-Buscalioni and P. Coveney 




10 20 30 40 50 60 " 
time [t] 



3L 




time [t] 



FIG.4. R. Delgado-Buscalioni and P. Coveney 




FIG. 5. R. Delgado-Buscalioni and P. Coveney 



0.5 0.6 0.7 0.8 0.9 1 

P [cr 3 ] 



FIG. 6. R. Delgado-Buscalioni and P. Coveney 




FIG. 7. R. Delgado-Buscalioni and P. Coveney 




FIG. 8. R. Delgado-Buscalioni and P. Coveney 



